

####################################################
####################################################
##########     Subsidiary functions       ##########
####################################################
####################################################

# The purpose of this file is to kee a log of relevant functions for my project.
# Headers of each section describe which file the functions are relevant for

 
## Extract the year from the quarter value
convert.YearQ.to.Year <- function(x) {
  y <- as.character(x)
  return(as.integer(paste0(substr(x = y, start = 1, stop = 4))))
} # End extracting the year

#---------------------------------------------------------------------------------------------

## Input normal color, output specified alpha-version
makeTransparent = function(..., alpha=0.5) {
  
  if(alpha<0 | alpha>1) stop("alpha must be between 0 and 1")
  
  alpha = floor(255*alpha)  
  newColor = col2rgb(col=unlist(list(...)), alpha=FALSE)
  
  .makeTransparent = function(col, alpha) {
    rgb(red=col[1], green=col[2], blue=col[3], alpha=alpha, maxColorValue=255)
  }
  
  newColor = apply(newColor, 2, .makeTransparent, alpha=alpha)
  
  return(newColor)
  
} # End function generating alpha-version of plots

#---------------------------------------------------------------------------------------------


## Extract slopes of functions outputted from the mixed-effect models
extract.slopes <- function(input.model, col.intercept, col.slope) {
  # Extract the main data-frame of interest
  input.df <- coef(input.model)$GEO
  # Generate a list of countries
  list.countries <- row.names(input.df)
  # Construct main data.frame to return
  return.df <- data.frame(countries = list.countries,
                          intercept.value = numeric(length = length(list.countries)),
                          slope.value = numeric(length= length(list.countries)))
  # Extract intercept values of countries
  return.df$intercept.value <- input.df[ , colnames(input.df) == col.intercept]
  # Extract slope values of countries
  return.df$slope.value <- input.df[ , colnames(input.df) == col.slope]
  # Return main data.frame
  return(return.df)
} # End function that extracts all the fixed- and marginal-coefficients

#---------------------------------------------------------------------------------------------


#### Scale data so I can interpret resulting beta coefficients
scale.data.frame <- function(input.data) {
  
  # Remove columns that aren't relevant
  do.not.keep <- c("Human_status", "Geneva", "Total_positive", "Rejected")
  input.data <- input.data[ , !(colnames(input.data) %in% do.not.keep)]
  # Define variables I'll always keep when partitioning
  keep.always <- c("TIME", "GEO", "CITIZEN")
  
  #-------------------------------------------------------------------------------------------
  # :::: CATEGORY 1: MEAN-VALUE FOR (GEO, CITIZEN) ::::
  # Keep only the 'received' category, which is unique for this type of normalization
  keep.1 <- c("Received")
  mean.geo.citizen <- input.data[ , colnames(input.data) %in% c(keep.always, keep.1) ]
  # Aggregate based on country
  attach(mean.geo.citizen)
  agg.mean.geo.citizen <- aggregate(x = mean.geo.citizen$Received, 
                                    by = list(GEO, CITIZEN), 
                                    FUN = mean, na.rm = TRUE)
  detach(mean.geo.citizen)
  colnames(agg.mean.geo.citizen) <- c("GEO", "CITIZEN", "Received.mean.raw")
  # Now loop through each citizen category and normalize accordingly
  agg.mean.geo.citizen$Received.mean.scale <- scale(x = agg.mean.geo.citizen$Received.mean.raw, center = TRUE, scale = TRUE)[,1]
  
  #---------# This code normalizes the applications received << per country >>
  agg.mean.geo.citizen$Received.mean.scale <- numeric(length = length(agg.mean.geo.citizen$GEO))
  for (c in unique(agg.mean.geo.citizen$CITIZEN)) {
    # Subset the vector to only include values from that CITIZEN group
    vec <- agg.mean.geo.citizen$Received.mean.raw[ agg.mean.geo.citizen$CITIZEN == c]
    agg.mean.geo.citizen$Received.mean.scale[ agg.mean.geo.citizen$CITIZEN == c] <- scale(x = vec, center = TRUE, scale = TRUE)[ , 1]
  } # End looping through each citizen category
  #---------# This code normalizes the applications received << per country >>
  
  ## Add to existing data.frame
  input.data <- merge(x = input.data, y = agg.mean.geo.citizen, 
                      by = c("GEO", "CITIZEN"), all.x = TRUE)
  # Remove unessesary items
  rm(agg.mean.geo.citizen,mean.geo.citizen,keep.1)
  #-------------------------------------------------------------------------------------------
  
  #-------------------------------------------------------------------------------------------
  # :::: CATEGORY 2: MEAN-VALUE FOR (GEO) ::::
  ## First, generate the mean-value vector, and normalize accordingly
  keep.2 <- c("NumArticles.Baseline", "NumArticles.Refugee", "NumArticles.Ratio", "MEDIA", "Tone.Baseline" , "Tone.Refugee",
              "Tone.Rel", "NumArticles.Rel", "CPI", "GDP", "Gov", "Unemploy", "Ideology", "Press.Freedom")
  # .. keep only one citizen category
  c <- "Turkey"
  time.invariant <- input.data[ input.data$CITIZEN == c , colnames(input.data) %in% c(keep.always, keep.2) ]
  # Next, aggregate the vector according to European country code 
  attach(time.invariant)
  small.step <- time.invariant[ , colnames(time.invariant) %in% keep.2 ]
  agg.time.invariant <- aggregate(x = small.step, # Keep numeric-relevant columns
                                  by = list(GEO), # aggregate only by the country
                                  FUN = mean,
                                  na.rm = TRUE) # and calculate the mean
  detach(time.invariant)
  ## Replace names appropriately!!! 
  #... first item is always GEO
  colnames(agg.time.invariant)[1] = "GEO"
  scale.time.invariant <- agg.time.invariant #... copy item before names are changed
  # ... now fix the rest
  colnames(agg.time.invariant) <- append.colname(names = colnames(agg.time.invariant), append.item = ".mean.raw")
  # Finally, normalize each column
  colnames(scale.time.invariant) <- append.colname(names = colnames(scale.time.invariant), append.item = ".mean.scale")
  for (col in keep.2) {
    vec <- agg.time.invariant[ , colnames(agg.time.invariant) == paste0(col, ".mean.raw") ]
    if (col %in% c("NumArticles.Ratio", "Tone.Baseline", "Tone.Refugee")) {
      scale.time.invariant[ , colnames(scale.time.invariant) == paste0(col, ".mean.scale")] = vec
    } else {
      scale.time.invariant[ , colnames(scale.time.invariant) == paste0(col, ".mean.scale")] = scale(x = vec, center = TRUE, scale = TRUE)[ , 1]
    }
  } # End looping through each column
  # Finally, merge data together
  merge.time.invariant <- merge(x = agg.time.invariant, y = scale.time.invariant, by = c("GEO"))
  input.data <- merge(x = input.data, y = merge.time.invariant, by = "GEO", all.x = TRUE)
  # Remove unecessary fluff from data storage
  rm(keep.2,c,col,vec,time.invariant,agg.time.invariant,scale.time.invariant,merge.time.invariant)
  #-------------------------------------------------------------------------------------------
  
  
  #-------------------------------------------------------------------------------------------
  # :::: CATEGORY 3: MEAN-VALUE FOR WITHIN-EFFECTS ::::
  keep.3 <- c("NumArticles.Baseline", "NumArticles.Refugee", "NumArticles.Ratio", "MEDIA", "Tone.Baseline" , "Tone.Refugee",
              "Tone.Rel", "NumArticles.Rel", "Received", "CPI", "GDP", "Gov", "Unemploy", "Ideology", "Press.Freedom")
  # List of column names for scaled vectors
  replace.keep.3 <- sapply(X = keep.3, FUN = function(x) {return(paste0(x, ".scale"))})
  # Generate new data-frame with scaled values
  scale.input.data <- input.data[ , colnames(input.data) %in% c(keep.always, keep.3) ]
  replace
  
  colnames(scale.input.data)[] <- c(colnames(scale.input.data)[1:3], replace.keep.3)
  
  # Loop through every country, citizen, and column to normalize everything appropriately
  for (country in unique(input.data$GEO)) {
    for (citizen in unique(input.data$CITIZEN)) {
      for (i in 1:length(keep.3)) {
        # Subset the global dataframe according to country and country-of-origin
        vec <- input.data[ input.data$GEO == country &
                             input.data$CITIZEN == citizen , colnames(input.data) == keep.3[i] ]
        # Next, replace the column by normalizing according to mean=0 and sd=1
        if (keep.3[i] %in% c("NumArticles.Ratio", "Tone.Baseline", "Tone.Refugee")) {
          scale.input.data[ scale.input.data$GEO == country & # match according to countr
                              scale.input.data$CITIZEN == citizen , # match according to citizen
                            colnames(scale.input.data) == replace.keep.3[i][[1]] ] = vec
        } else {
          scale.input.data[ scale.input.data$GEO == country & # match according to countr
                              scale.input.data$CITIZEN == citizen , # match according to citizen
                            colnames(scale.input.data) == replace.keep.3[i][[1]] ] = scale(x = vec, center = TRUE, scale = TRUE)[,1]
        }
      } # End looping through each column vector to normalize
    } # End looping through each 'country-of-origin'
  } # End looping through each country
  # Lastly, merge two vectors together
  input.data <- merge(x = input.data, y = scale.input.data, by = c("GEO", "CITIZEN", "TIME"), all.x = TRUE)
  # Remove superfluous data
  rm(citizen, country, i, keep.3, replace.keep.3, vec, scale.input.data)
  #-------------------------------------------------------------------------------------------
  
  ## Lastly, return the cleaned & normalized data structure
  return(input.data)
} # End scaling data --- to acquire regression coefficients

#---------------------------------------------------------------------------------------------


## -->> Function for appending strings to end of column name
append.colname <- function(names, append.item) {
  # Loop through every item in the column.names
  for (i in 1:length(names)) {
    if (!(names[i] %in% c("GEO", "TIME", "CITIZEN"))) {
      names[i] = paste0(names[i], append.item)
    } # End checking if the value is not a categorical variable
  } # End looping through every column name
  return(names)
} # End function that appends items to end of column names in dataframe

#---------------------------------------------------------------------------------------------

## Function specifically to clean TIME data of CPI
clean.CPI <- function(x) {
  x.return <- paste0(substr(x, 1, 4),
                     substr(x, 6, 7))
  return(x.return)
}

#---------------------------------------------------------------------------------------------

## Quick function to quickly clean my data format
clean.economic.files <- function(imported.file) {
  require(zoo)
  # Convert ":" character into NA -- and change from string to integers
  imported.file$Value <- as.character(imported.file$Value)
  for (i in 1:length(imported.file$Value)) {
    imported.file$Value[i] = gsub(pattern = ",", replacement = "", x = imported.file$Value[i])
  }
  imported.file[imported.file == ":"] <- NA
  imported.file$Value <- as.numeric(imported.file$Value)
  
  # Finally, convert "GEO" and "CITIZEN" to the character class
  imported.file$GEO <- as.character(imported.file$GEO)
  imported.file[imported.file == "Germany (until 1990 former territory of the FRG)"] <- "Germany"
  
  # Convert "TIME" into an actual "data"-compatible format
  imported.file$TIME <- sapply(X = as.character(imported.file$TIME), FUN = add.space.TIME)
  build.vector <- c()
  for (t in imported.file$TIME) {
    build.vector <- c( as.yearqtr(t) , build.vector)
  }
  imported.file$TIME <- rev(build.vector)
  
  # Lastly, clean the imported file
  return(imported.file)
}

#---------------------------------------------------------------------------------------------


## Quick function to quickly clean my data format
clean.received.files <- function(imported.file) {
  # Convert ":" character into NA -- and change from string to integers
  imported.file$Value <- as.character(imported.file$Value)
  imported.file[imported.file == ":"] <- "0"
  imported.file$Value <- as.integer(gsub(",", "", imported.file$Value))
  imported.file$Value[is.na(imported.file$Value)] = 0
  
  # Finally, convert "GEO" and "CITIZEN" to the character class
  imported.file$GEO <- as.character(imported.file$GEO)
  imported.file$CITIZEN <- as.character(imported.file$CITIZEN)
  imported.file[imported.file == "Germany (until 1990 former territory of the FRG)"] <- "Germany"
  imported.file[imported.file == "China (including Hong Kong)"] <- "China"
  imported.file[imported.file == "Kosovo (under United Nations Security Council Resolution 1244/99)"] <- "Kosovo"
  imported.file[imported.file == "Democratic Republic of the Congo"] <- "DRC"
  
  # Generate a dummy column specifying these applications were received by the government
  imported.file$DECISION <- rep(x = "Received", times = length(imported.file$GEO))
  
  # Convert "TIME" into an actual "data"-compatible format
  imported.file$TIME <- as.character(imported.file$TIME)
  imported.file$TIME <- convert.MONTH.to.YearQ(TIME.var = imported.file$TIME)
  
  # Lastly, clean the imported file
  return(imported.file)
}

#---------------------------------------------------------------------------------------------


## Quick function to quickly clean my data format
clean.decision.files <- function(imported.file) {
  # Convert ":" character into NA -- and change from string to integers
  imported.file$Value <- as.character(imported.file$Value)
  imported.file[imported.file == ":"] <- "0"
  imported.file$Value <- as.integer(gsub(",", "", imported.file$Value))
  imported.file$Value[is.na(imported.file$Value)] = 0
  
  # Finally, convert "GEO" and "CITIZEN" to the character class
  imported.file$GEO <- as.character(imported.file$GEO)
  imported.file$DECISION <- as.character(imported.file$DECISION)
  imported.file$CITIZEN <- as.character(imported.file$CITIZEN)
  imported.file[imported.file == "Germany (until 1990 former territory of the FRG)"] <- "Germany"
  imported.file[imported.file == "China (including Hong Kong)"] <- "China"
  imported.file[imported.file == "Kosovo (under United Nations Security Council Resolution 1244/99)"] <- "Kosovo"
  imported.file[imported.file == "Democratic Republic of the Congo"] <- "DRC"
  
  # Convert "TIME" into an actual "data"-compatible format
  imported.file$TIME <- as.character(imported.file$TIME)
  if ( substr(x = imported.file$TIME[1], start = 5, stop = 5 ) == "M" ) {
    imported.file$TIME <- convert.MONTH.to.YearQ(TIME.var = imported.file$TIME)
  } else if ( substr(x = imported.file$TIME[1], start = 5, stop = 5 ) == "Q" ) {
    imported.file$TIME <- convert.QUARTER.to.YearQ(TIME.var = imported.file$TIME)
  }
  
  # Lastly, clean the imported file
  return(imported.file)
}

#---------------------------------------------------------------------------------------------


## Quick function to convert EU data to appropriate "date" formate
convert.MONTH.to.YearQ <- function(TIME.var) {
  require(zoo)
  return.TIME <- as.Date( paste0( substr(TIME.var, 1 , 4), "-",
                                  substr(TIME.var, 6 , 7), "-01") )
  return(as.yearqtr(return.TIME))
}

#---------------------------------------------------------------------------------------------

## Quick function to convert EU data to appropriate "date" formate
convert.QUARTER.to.YearQ <- function(TIME.var) {
  require(zoo)
  return(as.yearqtr(TIME.var))
}

#---------------------------------------------------------------------------------------------


## Add a space between year and quarter time
add.space.TIME <- function(TIME.var) {
  return.TIME <- paste0(substr(x = TIME.var, start = 1, stop = 4),
                        " ",
                        substr(x = TIME.var, start = 5, stop = 6))
  return(return.TIME)
}

#---------------------------------------------------------------------------------------------

## Quick function to convert EU data to appropriate "date" formate
aggregate.asylum.apps <- function(input) {
  # Aggregate the vector with function "sum"
  attach(input)
  return.val <- aggregate(x = input$Value, 
                          by = list(TIME, GEO, CITIZEN, DECISION), 
                          FUN = sum, 
                          na.omit = TRUE)
  detach(input)
  colnames(return.val) <- c("TIME", "GEO", "CITIZEN", "DECISION", "Value")
  # Lastly, return the constructed data-frame
  return(return.val)
} # End function for aggregating monthly to quarterly data


#---------------------------------------------------------------------------------------------


merge.asylum.data.frames <- function(asylum.apps) {
  
  received.df <- asylum.apps[asylum.apps$DECISION == "Received" , ][ , c(1:3,5)]
  colnames(received.df) <- c("TIME", "GEO", "CITIZEN", "Received")
  
  tp.df <- asylum.apps[asylum.apps$DECISION == "Total positive decisions" , ][ , c(1:3,5)]
  colnames(tp.df) <- c("TIME", "GEO", "CITIZEN", "Total_positive")
  
  human.df <- asylum.apps[asylum.apps$DECISION == "Humanitarian status" , ][ , c(1:3,5)]
  colnames(human.df) <- c("TIME", "GEO", "CITIZEN", "Human_status")
  
  geneva.df <- asylum.apps[asylum.apps$DECISION == "Geneva Convention status" , ][ , c(1:3,5)]
  colnames(geneva.df) <- c("TIME", "GEO", "CITIZEN", "Geneva")
  
  rejected.df <- asylum.apps[asylum.apps$DECISION == "Rejected" , ][ , c(1:3,5)]
  colnames(rejected.df) <- c("TIME", "GEO", "CITIZEN", "Rejected")
  
  # Put sub-dataframes together in a list  
  list.apps <- list(received.df, tp.df, human.df, geneva.df, rejected.df)
  # Now merge all the lists together
  merged.asylum.apps <- Reduce(function(...) merge(..., by = c("TIME", "GEO", "CITIZEN"), all=T), list.apps)
  
  # Calculate the acceptance rate (total positive) / (total positive + rejected)
  merged.asylum.apps$Acceptance.Rate <- merged.asylum.apps$Total_positive / (merged.asylum.apps$Total_positive + merged.asylum.apps$Rejected)
  
  return(merged.asylum.apps)
  
} # End merging all of the data frames together

#---------------------------------------------------------------------------------------------

### Interchange the country name with the country code
exchange.CountryName.to.Code <- function(name.vector, desired.output) {
  country.code <- c("FR", "SW", "UK", "SZ", "AU", "BE", "BU", "HR", "EZ", "DA", "EN", "FI", 
                    "GR", "HU", "IT", "LG", "LS", "LH", "NL", "NO", "PL", "PO",
                    "RO", "LO", "SP", "GM")
  country.name <- c("France", "Sweden", "United Kingdom", "Switzerland", "Austria", 
                    "Belgium", "Bulgaria", "Croatia", "Czech Republic", "Denmark", "Estonia",
                    "Finland", "Greece", "Hungary", "Italy", "Latvia", "Liechtenstein",
                    "Lithuania", "Netherlands", "Norway", "Poland", "Portugal", "Romania", 
                    "Slovakia", "Spain", "Germany")
  # Loop through every "name" in vector and exchange for appropriate item
  if (desired.output == "CountryName") {
    for (i in 1:length(name.vector)) {
      name.vector[i] = country.name[ which(name.vector[i] == country.code) ]
    } # End looping through every country code to be exchanged
  } else if (desired.output == "CountryCode") {
    for (i in 1:length(name.vector)) {
      name.vector[i] = country.code[ which(name.vector[i] == country.name) ]
    } # End looping through every country name to be exchanged
  } # End determining of desired output is "CountryName" or "CountryCode"
  return(name.vector)
} # End exchanging country code for country name

#---------------------------------------------------------------------------------------------


## Quick function to quickly clean my data format
clean.gdelt.file <- function(imported.file) {  
  require(zoo)
  # Change column names
  colnames(imported.file) <- c("TIME", "Tone.Baseline", "NumArticles.Baseline",
                               "GEO", "Tone.Refugee", "NumArticles.Refugee",
                               "Tone.Rel", "NumArticles.Rel")
  # Convert country characters from "factor" to "character"
  imported.file$GEO <- as.character(imported.file$GEO)
  
  # Convert "TIME" into an actual "data"-compatible format
  imported.file$TIME <- as.character(imported.file$TIME)
  build.vector <- c()
  for (t in imported.file$TIME) {
    build.vector <- c( as.yearqtr(as.Date(t)) , build.vector)
  }
  imported.file$TIME <- rev(build.vector)
  
  # Generate a media file
  imported.file$MEDIA <- imported.file$Tone.Rel * imported.file$NumArticles.Rel
  
  ## Lastly, aggregate data
  # --> first aggregate according to variables that need to be summed
  agg.sum <- aggregate(cbind(NumArticles.Baseline, NumArticles.Refugee, MEDIA) ~ TIME + GEO,
                       data = imported.file, FUN = sum, na.rm = TRUE)
  # --> second aggregate according to variables via mean
  agg.mean <- aggregate(cbind(Tone.Baseline, Tone.Refugee, Tone.Rel, NumArticles.Rel) ~ TIME + GEO,
                        data = imported.file, FUN = mean, na.rm = TRUE)
  agg.mean
  # --> third merge the two dataframes together
  imported.file <- merge(agg.sum, agg.mean, by = c("GEO", "TIME"))
  
  # Lastly, clean the imported file
  return(imported.file)
}

#---------------------------------------------------------------------------------------------


## Quick function to quickly clean my data format
clean.files.monthly <- function(imported.file, input.decision.type) {
  # Convert ":" character into NA -- and change from string to integers
  imported.file$Value <- as.character(imported.file$Value)
  imported.file[imported.file == ":"] <- NA
  imported.file$Value <- as.integer(gsub(",", "", imported.file$Value))
  
  # Finally, convert "GEO" and "CITIZEN" to the character class
  imported.file$GEO <- as.character(imported.file$GEO)
  imported.file$CITIZEN <- as.character(imported.file$CITIZEN)
  imported.file[imported.file == "Germany (until 1990 former territory of the FRG)"] <- "Germany"
  imported.file[imported.file == "China (including Hong Kong)"] <- "China"
  imported.file[imported.file == "Kosovo (under United Nations Security Council Resolution 1244/99)"] <- "Kosovo"
  imported.file[imported.file == "Democratic Republic of the Congo"] <- "DRC"
  
  # Generate a dummy column specifying these applications were received by the government
  colnames(imported.file)[ colnames(imported.file) == "Value" ] = input.decision.type
  #imported.file$DECISION <- rep(x = input.decision.type, times = length(imported.file$GEO))
  
  # Convert "TIME" into an actual "data"-compatible format
  imported.file$TIME <- as.character(imported.file$TIME)
  imported.file$TIME <- convert.string.to.Date(value = imported.file$TIME) 
  
  # Lastly, clean the imported file
  return(imported.file)
} # End cleaning imported-file for montly data

#---------------------------------------------------------------------------------------------

## Quick function to quickly clean my data format
clean.gdelt.file.montly <- function(imported.file) {  
  require(zoo)
  # Change column names
  colnames(imported.file) <- c("TIME", "Tone.Baseline", "NumArticles.Baseline",
                               "GEO", "Tone.Refugee", "NumArticles.Refugee",
                               "Tone.Rel", "NumArticles.Rel")
  # Convert country characters from "factor" to "character"
  imported.file$GEO <- as.character(imported.file$GEO)
  
  # Convert "TIME" into an actual "data"-compatible format
  imported.file$TIME <- convert.string.to.Date(value = as.character(imported.file$TIME))
  
  # Generate a media file
  imported.file$MEDIA <- imported.file$Tone.Rel * imported.file$NumArticles.Rel
  
  # Lastly, clean the imported file
  return(imported.file)
} # End function cleaning all the GDELT files

#---------------------------------------------------------------------------------------------


### Subordinate function to "clean.received.files.monthly"
convert.string.to.Date <- function(value) {
  x <- paste0(substr(x = value, start = 1, stop = 4),
              "-",
              substr(x = value, start = 6, stop = 7),
              "-01")
  return(as.Date(x))
} # End function converting all string characters to dates

#---------------------------------------------------------------------------------------------


### Interchange the country name with the country code
exchange.CountryName.to.Code <- function(name.vector, desired.output) {
  country.code <- c("FR", "SW", "UK", "SZ", "AU", "BE", "BU", "HR", "EZ", "DA", "EN", "FI", 
                    "GR", "HU", "IT", "LG", "LS", "LH", "NL", "NO", "PL", "PO",
                    "RO", "LO", "SP", "GM")
  country.name <- c("France", "Sweden", "United Kingdom", "Switzerland", "Austria", 
                    "Belgium", "Bulgaria", "Croatia", "Czech Republic", "Denmark", "Estonia",
                    "Finland", "Greece", "Hungary", "Italy", "Latvia", "Liechtenstein",
                    "Lithuania", "Netherlands", "Norway", "Poland", "Portugal", "Romania", 
                    "Slovakia", "Spain", "Germany")
  # Loop through every "name" in vector and exchange for appropriate item
  if (desired.output == "CountryName") {
    for (i in 1:length(name.vector)) {
      name.vector[i] = country.name[ which(name.vector[i] == country.code) ]
    } # End looping through every country code to be exchanged
  } else if (desired.output == "CountryCode") {
    for (i in 1:length(name.vector)) {
      name.vector[i] = country.code[ which(name.vector[i] == country.name) ]
    } # End looping through every country name to be exchanged
  } # End determining of desired output is "CountryName" or "CountryCode"
  return(name.vector)
} # End exchanging country code for country name

#---------------------------------------------------------------------------------------------


### Function that calculates ranking for country-of-origin
ranking.df <- function(final.data, max.rank = 30) {
  
  description.data <- final.data[ !(final.data$CITIZEN %in% c("Extra EU-28", "European Economic Area (EEA18-2004, EEA28-2006, EEA30-2013, EEA31)")) , ]
  description.data$YEAR <- convert.YearQ.to.Year(x = description.data$TIME)
  #####################################################################################
  #     Determine the total number of asylum decisions made per country 
  #####################################################################################
  
  return.df <- data.frame(rank = 1:max.rank) 
  
  #-------------------------------------------------------------------------------------------
  ######### 2002 -- 2007 ::: and keep only "total positive" & "rejected"
  data.2002.2007 <- description.data[ description.data$YEAR > 2001 & description.data$YEAR < 2008 , ]
  # Aggregate data
  attach(data.2002.2007)
  agg.2002.2007 <- aggregate(x = cbind(data.2002.2007$Total_positive, data.2002.2007$Rejected),
                             by = list(GEO, CITIZEN),
                             FUN = sum)
  detach(data.2002.2007)
  colnames(agg.2002.2007) <- c("GEO", "CITIZEN", "Total_positive", "Rejected")
  agg.2002.2007$Total <- agg.2002.2007$Total_positive + agg.2002.2007$Rejected
  # Determine which "country of origin"s are most relevant
  attach(agg.2002.2007)
  popular.2002.2007 <- aggregate(x = agg.2002.2007$Total, 
                                 by = list(CITIZEN), 
                                 FUN = sum)
  detach(agg.2002.2007)
  popular.2002.2007 <- popular.2002.2007[ order(popular.2002.2007$x, decreasing = TRUE) , ]  
  ranking.2002.2007 <- popular.2002.2007$Group.1
  rm(data.2002.2007)
  return.df$Ranking.02.07 <- ranking.2002.2007[1:max.rank]
  #-------------------------------------------------------------------------------------------
  
  
  #-------------------------------------------------------------------------------------------
  ######### 2008 -- 2016 ::: and keep only "total positive" & "rejected"
  data.2008.2016 <- description.data[ description.data$YEAR > 2007 & description.data$YEAR < 2017 , ]
  # Aggregate data
  attach(data.2008.2016)
  agg.2008.2016 <- aggregate(x = cbind(data.2008.2016$Total_positive, data.2008.2016$Rejected),
                             by = list(GEO, CITIZEN),
                             FUN = sum)
  detach(data.2008.2016)
  colnames(agg.2008.2016) <- c("GEO", "CITIZEN", "Total_positive", "Rejected")
  agg.2008.2016$Total <- agg.2008.2016$Total_positive + agg.2008.2016$Rejected
  # Determine which "country of origin"s are most relevant
  attach(agg.2008.2016)
  popular.2008.2016 <- aggregate(x = agg.2008.2016$Total, 
                                 by = list(CITIZEN), 
                                 FUN = sum)
  detach(agg.2008.2016)
  popular.2008.2016 <- popular.2008.2016[ order(popular.2008.2016$x, decreasing = TRUE) , ]
  ranking.2008.2016 <- popular.2008.2016$Group.1
  rm(data.2008.2016)
  return.df$Ranking.08.16 <- ranking.2008.2016[1:max.rank]
  #-------------------------------------------------------------------------------------------
  
  
  #-------------------------------------------------------------------------------------------
  ######### 2002 -- 2016 ::: and keep only "total positive" & "rejected"
  data.2002.2016 <- description.data[ description.data$YEAR > 2001 & description.data$YEAR < 2017 , ]
  # Aggregate data
  attach(data.2002.2016)
  agg.2002.2016 <- aggregate(x = cbind(data.2002.2016$Total_positive, data.2002.2016$Rejected),
                             by = list(GEO, CITIZEN),
                             FUN = sum)
  detach(data.2002.2016)
  colnames(agg.2002.2016) <- c("GEO", "CITIZEN", "Total_positive", "Rejected")
  agg.2002.2016$Total <- agg.2002.2016$Total_positive + agg.2002.2016$Rejected
  # Determine which "country of origin"s are most relevant
  attach(agg.2002.2016)
  popular.2002.2016 <- aggregate(x = agg.2002.2016$Total, 
                                 by = list(CITIZEN), 
                                 FUN = sum)
  detach(agg.2002.2016)
  popular.2002.2016 <- popular.2002.2016[ order(popular.2002.2016$x, decreasing = TRUE) , ]
  ranking.2002.2016 <- popular.2002.2016$Group.1
  rm(data.2002.2016)
  return.df$Ranking.02.16 <- ranking.2002.2016[1:max.rank]
  #-------------------------------------------------------------------------------------------
  return(return.df)
} ## End function that computes rankings for country-of-origin (across Europe)

### Function that re-organizes data for bean-plot
reorganize.beanplot.data <- function(plot.data, country.order.CODE) {
  
  for (c in country.order.CODE) {
    if ( !(is.na(c)) ) {
      if (c == country.order.CODE[1]) {
        return.data.frame <- rbind(plot.data[ plot.data$CITIZEN == c & plot.data$Section == "2002-07" ,  ],
                                   plot.data[ plot.data$CITIZEN == c & plot.data$Section == "2008-16" ,  ])
      } else {
        return.data.frame <- rbind(return.data.frame,
                                   plot.data[ plot.data$CITIZEN == c & plot.data$Section == "2002-07" ,  ],
                                   plot.data[ plot.data$CITIZEN == c & plot.data$Section == "2008-16" ,  ])
      } # End looping through order of countries
    }  # End checking if element is NA --- not a country!
  } # End looping through every unique country in list
  
  return.data <- return.data.frame
  return.data$CITIZEN <- as.character(return.data$CITIZEN)
  return.data$Section <- as.character(return.data$Section)
  for (i in 1:length(return.data$Section)) {
    if (return.data$Section[i] == "2002-07") {
      return.data$CITIZEN[i] = paste0(return.data$CITIZEN[i], " 1")
    } else {
      return.data$CITIZEN[i] = paste0(return.data$CITIZEN[i], " 2")
    } # Otherwise continue
  } # End looping through entire section
  
  return(return.data)
}

#---------------------------------------------------------------------------------------------


p.num <- function(x, digits) {
  x.val <- x[[1]]
  x.p <- x[[5]]
  
  x.round <- as.character(round(x = x.val, digits = digits))
  if (x.p > 0.01 & x.p < 0.05) {
    x.round <- paste0(x.round,"^{*}")
  } else if (x.p > 0.001 & x.p < 0.01) {
    x.round <- paste0(x.round,"^{**}")
  } else if (x.p < 0.001) {
    x.round <- paste0(x.round,"^{***}")
  } #... end checking if coefficient is significant
  return(x.round)
}

#---------------------------------------------------------------------------------------------



## Small function calculating interesting model results
print.model.results <- function(input.models, CONTROLS = FALSE, MEDIA = FALSE) {
  
  # Print main summary values of coefficients!!!
  print(testEstimates(input.models, var.comp = TRUE))
  
  ### ---- General model statistics
  # Print the number of observations
  print(paste0("Number of observations: ", input.models$imp1@devcomp$dims[1]))
  # Calculate R2 values
  r.squared.results <- as.data.frame(t(sapply(X = input.models, FUN =r.squaredGLMM)))
  # Print Fixed-effect R^2 results
  print(paste0("Fixed-effect R^2: ", poolR2(input.R = r.squared.results[ , 1 ]))) #--- fixed-effect R^2
  # Print conditional R^2 results
  print(paste0("Conditional R^2: ", poolR2(input.R= r.squared.results[ , 2 ]))) #--- conditional R^2
  # Print AICc result
  print(paste0("AICc: ", mean(unlist(lapply(input.models, FUN = AICc))))) #--- AICc value, corrected AIC estimation now accounting for finite sample size
  print(paste0("BIC: ", mean(unlist(lapply(input.models, FUN = AICc))))) #--- AICc value, corrected AIC estimation now accounting for finite sample size
  
  corr.df <- data.frame(intercept.NumArticles = numeric(length = length(input.models)),
                        intercept.Tone = numeric(length = length(input.models)),
                        NumArticles.Tone = numeric(length = length(input.models)),
                        intercept.Received.norm = numeric(length = length(input.models)))
  ### Calculate correlative values!! 
  ## Correlations between coefficients
  require(parallel)
  summary.list <- mclapply(input.models, FUN = summary, mc.cores = 2)
  for (i in 1:length(input.models)) {
    test <- summary.list[[i]]
    
    ### --- >> Media variables
    if (MEDIA == TRUE) {
      corr.df$intercept.NumArticles[i] = attributes(test$varcor$GEO)$correlation[ 1 , 2]
      corr.df$intercept.Tone[i] = attributes(test$varcor$GEO)$correlation[1 , 3]
      corr.df$NumArticles.Tone[i] = attributes(test$varcor$GEO)$correlation[2 , 3]
    } # End checking if media variables exist
    
    ### --- >> Control variables
    if (CONTROLS == TRUE) {
      corr.df$intercept.Received.norm[i] = attributes(test$varcor$CITIZEN.TIME)$correlation[ 1 , 2]
    } # End checking if control variables exist
    
  } # End looping through each model
  
  ### ---- >>> PRINT RESULTS
  if ( MEDIA == TRUE) {
    print("::: MEDIA :::")
    print(paste0("Cor(Intercept, NumArticles): ", poolR2(input.R = corr.df$intercept.NumArticles)))
    print(paste0("Cor(Intercept, Tone): ", poolR2(input.R = corr.df$intercept.Tone)))
    print(paste0("Cor(NumArticles, Tone): ", poolR2(input.R = corr.df$NumArticles.Tone)))
  } # End checking if media variables exist
  
  ## Correlations between coefficients
  if ( CONTROLS == TRUE ) {
    print("::: CITIZEN.TIME :::") 
    print(paste0("Cor(Intercept, Received.norm): ", poolR2(input.R = corr.df$intercept.Received.norm)))
  } ## End checking if model has correlative values
  
} # End printing results



#---------------------------------------------------------------------------------------------


#########################################################################
###  Write 'heirarchical regression' function
#########################################################################

### --->>> Function for pooling R-squared values
poolR2 <- function(input.R) {
  z.value <- unlist(fisherz(input.R))
  mean.z <- mean(z.value)
  return( fisherz2r(mean.z) )
} ## End function pooling R-squared values

#---------------------------------------------------------------------------------------------

clean.data <- function(input.data) {
  ## Conduct data cleaning
  input.data$Acceptance.Rate[ input.data$Acceptance.Rate == 1.0 ] = 0.999
  input.data$Acceptance.Rate[ input.data$Acceptance.Rate == 0.0 ] = 0.001
  input.data$Logit.NumArticles <- log(input.data$NumArticles.Ratio.scale / (1-input.data$NumArticles.Ratio.scale))
  input.data$Logit.NumArticles.mean <- log(input.data$NumArticles.Ratio.mean.scale / (1-input.data$NumArticles.Ratio.mean.scale))
  
  input.data$Logit.Tone.R <- log((100 - input.data$Tone.Refugee.scale) / (input.data$Tone.Refugee.scale + 100))
  input.data$Logit.Tone.B <- log((100 - input.data$Tone.Baseline.scale) / (input.data$Tone.Baseline.scale + 100))
  input.data$Logit.Tone.mean.R <- log((100 - input.data$Tone.Refugee.mean.scale) / (input.data$Tone.Refugee.mean.scale + 100))
  input.data$Logit.Tone.mean.B <- log((100 - input.data$Tone.Baseline.mean.scale) / (input.data$Tone.Baseline.mean.scale + 100))
  
  #... define dependent variable
  input.data$Logit.Acceptance <- log((input.data$Acceptance.Rate) / (1-input.data$Acceptance.Rate))
  #... define time as dummy variable
  input.data$TIME <- as.character(input.data$TIME)
  
  return(input.data)
}

#---------------------------------------------------------------------------------------------

## HALF REGRESSION --- MEDIA VARIABLES
lmer.regress.MEDIA <- function(input.data) {
  input.data <- clean.data(input.data = input.data)
  lmer.out <- lmer(formula = Acceptance.Rate ~ Logit.NumArticles + 
                     Logit.NumArticles.mean + 
                     Logit.Tone.R +
                     Logit.Tone.mean.R + 
                     (1 | CITIZEN.TIME) + (1 | GEO.CITIZEN),
                   data = input.data, REML = FALSE)
  #summary(lmer.out)
  #r.squaredGLMM(object = lmer.out)
  return(lmer.out)
} ## End lmer test formulation

#---------------------------------------------------------------------------------------------

## HALF REGRESSION --- MEDIA VARIABLES
lmer.regress.MEDIA.2 <- function(input.data) {
  input.data <- clean.data(input.data = input.data)
  lmer.out <- lmer(formula = Acceptance.Rate ~ Logit.NumArticles +
                     Logit.NumArticles.mean + 
                     I(Logit.Tone.R - Logit.Tone.B) + 
                     Logit.Tone.mean.R +
                     (1 | CITIZEN.TIME) + (1 | GEO.CITIZEN),
                   data = input.data, REML = FALSE)
  return(lmer.out)
} ## End lmer test formulation

#---------------------------------------------------------------------------------------------

## HALF REGRESSION --- RECEIVED VARIABLES
lmer.regress.RECEIVED <- function(input.data) {
  input.data <- clean.data(input.data = input.data)
  lmer.out <- lmer(formula = Logit.Acceptance ~ Received.scale + Received.mean.scale +
                     (1 | CITIZEN.TIME) + (1 | GEO.CITIZEN),
                   data = input.data, REML = FALSE)
  return(lmer.out)
} ## End lmer test formulation

#---------------------------------------------------------------------------------------------

## HALF REGRESSION --- CONTROL VARIABLES
lmer.regress.CONTROL <- function(input.data) {
  input.data <- clean.data(input.data = input.data)
  lmer.out <- lmer(formula = Logit.Acceptance ~ Received.scale + Received.mean.scale +
                     GDP.scale + GDP.mean.scale +
                     Gov.scale + Gov.mean.scale +
                     Unemploy.scale + Unemploy.mean.scale +
                     CPI.scale + CPI.mean.scale + 
                     Ideology.scale + Ideology.mean.scale +
                     Press.Freedom.scale + Press.Freedom.mean.scale +
                     (1 | CITIZEN.TIME) + (1 | GEO.CITIZEN),
                   data = input.data, REML = FALSE)
  return(lmer.out)
} ## End lmer test formulation

#---------------------------------------------------------------------------------------------

## FULL REGRESSION --- MEDIA + CONTROL VARIABLES
lmer.regress.FULL <- function(input.data) {
  input.data <- clean.data(input.data = input.data)
  lmer.out <- lmer(formula =  Logit.Acceptance ~ Logit.NumArticles + Logit.NumArticles.mean + 
                     Logit.Tone.R +
                     Logit.Tone.mean.R +
                     Received.scale + Received.mean.scale +
                     GDP.scale + GDP.mean.scale +
                     Gov.scale + Gov.mean.scale +
                     Unemploy.scale + Unemploy.mean.scale +
                     CPI.scale + CPI.mean.scale + 
                     Ideology.scale + Ideology.mean.scale +
                     Press.Freedom.scale + Press.Freedom.mean.scale +
                     (1 | CITIZEN.TIME) + (1 | GEO.CITIZEN),
                   data = input.data, REML = FALSE)
  return(lmer.out)
} ## End lmer test formulation

#---------------------------------------------------------------------------------------------

## FULL REGRESSION --- MEDIA + CONTROL VARIABLES
lmer.regress.FULL.2 <- function(input.data) {
  input.data <- clean.data(input.data = input.data)
  lmer.out <- lmer(formula =  Logit.Acceptance ~ Logit.NumArticles + Logit.NumArticles.mean + 
                     I(Logit.Tone.R - Logit.Tone.B) + 
                     Logit.Tone.mean.R +
                     Received.scale + Received.mean.scale +
                     GDP.scale + GDP.mean.scale +
                     Gov.scale + Gov.mean.scale +
                     Unemploy.scale + Unemploy.mean.scale +
                     CPI.scale + CPI.mean.scale + 
                     Ideology.scale + Ideology.mean.scale +
                     Press.Freedom.scale + Press.Freedom.mean.scale +
                     (1 | CITIZEN.TIME) + (1 | GEO.CITIZEN),
                   data = input.data, REML = FALSE)
  return(lmer.out)
} ## End lmer test formulation

#---------------------------------------------------------------------------------------------



out.of.sample.predictions <- function(num, input.full.data, model.type) {
  
  # Clean the input data in the same way as when the previous regressions were fitted
  input.clean.data <- clean.data(input.data = input.full.data)
  
  # Keep track of indices
  index.vec <- 1:(dim(input.clean.data)[1])
  
  # Sample for training data
  training.index <- sample(x = index.vec, size = round(x = length(index.vec)*0.8), replace = FALSE)
  
  # Now keep testing data (which is not in the training data)
  testing.index <- index.vec[ !(index.vec %in% training.index) ]
  
  
  #-------------# Fit model to the appropriate type
  if (model.type == 1) {
    
    model.fit <- lmer(formula = Acceptance.Rate ~ -1 +
                        Logit.NumArticles + 
                        Logit.NumArticles.mean + 
                        Logit.Tone.R + 
                        Logit.Tone.mean.R +
                        (1 | GEO) + (1 | CITIZEN) + (1 | TIME),
                      data = input.clean.data[training.index, ], REML = FALSE)
    
  } else if (model.type == 2) {
    
    model.fit <- lmer(formula = Acceptance.Rate ~ -1 +
                        Logit.NumArticles + Logit.NumArticles.mean + 
                        I(Logit.Tone.R - Logit.Tone.B) + 
                        Logit.Tone.mean.R +
                        (1 | GEO) + (1 | CITIZEN) + (1 | TIME),
                      data = input.clean.data[training.index, ], REML = FALSE)
    
  } else if (model.type == 3) {
    
    model.fit <- lmer(formula = Logit.Acceptance ~ -1 +
                        Received.scale + 
                        Received.mean.scale + 
                        (1 | GEO) + (1 | CITIZEN) + (1 | TIME),
                      data = input.clean.data[training.index, ], REML = FALSE)
    
    
  } else if (model.type == 4) {
    
    model.fit <- lmer(formula = Logit.Acceptance ~ -1 +
                        Received.scale + Received.mean.scale +
                        GDP.scale + GDP.mean.scale +
                        Gov.scale + Gov.mean.scale +
                        Unemploy.scale + Unemploy.mean.scale +
                        CPI.scale + CPI.mean.scale + 
                        Ideology.scale + Ideology.mean.scale +
                        Press.Freedom.scale + Press.Freedom.mean.scale +
                        (1 | GEO) + (1 | CITIZEN) + (1 | TIME),
                      data = input.clean.data[training.index, ], REML = FALSE)
    
  } else if (model.type == 5) {
    
    model.fit <- lmer(formula =  Logit.Acceptance ~ -1 +
                        Logit.NumArticles + Logit.NumArticles.mean + 
                        Logit.Tone.R +
                        Logit.Tone.mean.R +
                        Received.scale + Received.mean.scale +
                        GDP.scale + GDP.mean.scale +
                        Gov.scale + Gov.mean.scale +
                        Unemploy.scale + Unemploy.mean.scale +
                        CPI.scale + CPI.mean.scale + 
                        Ideology.scale + Ideology.mean.scale +
                        Press.Freedom.scale + Press.Freedom.mean.scale +
                        (1 | GEO) + (1 | CITIZEN) + (1 | TIME),
                      data = input.clean.data[training.index, ], REML = FALSE)
    
  } else if (model.type == 6) {
    
    model.fit <- lmer(formula =  Logit.Acceptance ~ -1 +
                        Logit.NumArticles + Logit.NumArticles.mean + 
                        I(Logit.Tone.R - Logit.Tone.B) + 
                        Logit.Tone.mean.R +
                        Received.scale + Received.mean.scale +
                        GDP.scale + GDP.mean.scale +
                        Gov.scale + Gov.mean.scale +
                        Unemploy.scale + Unemploy.mean.scale +
                        CPI.scale + CPI.mean.scale + 
                        Ideology.scale + Ideology.mean.scale +
                        Press.Freedom.scale + Press.Freedom.mean.scale +
                        (1 | GEO) + (1 | CITIZEN) + (1 | TIME),
                      data = input.clean.data[training.index, ], REML = FALSE)
  } else {
    model.fit <- NULL
  }  #----------------------------------------#
  
  #... keep old acceptance rates
  old.value <- input.clean.data[testing.index, ]$Acceptance.Rate
  old.value <- exp(old.value) / (exp(old.value) + 1)
  #... keep 
  predict.value <- predict(object = model.fit, input.clean.data[testing.index, ], allow.new.levels = TRUE)
  predict.value <- exp(predict.value) / (exp(predict.value) + 1)
  ### Utilize these values (have complete cases)
  input.values <- cbind(old.value, predict.value)
  #... the '100*' converts to percentage value
  input.values <- 100*input.values[complete.cases(input.values) , ]

  # Calculate Mean-square error (MSE)
  return.MSE <- MSE(y_pred = input.values[ , 2], y_true = input.values[ , 1])

  # Calculate Mean absolute value (MAV)
  return.MAV <- mean(abs(input.values[,2] - input.values[,1]))

  # Return results
  return(c(return.MSE,return.MAV))
} #... end "out.of.sample.predictions"

#---------------------------------------------------------------------------------------------






####################################################################################################################################
####################################################################################################################################
#    Main functions for granger-causality analysis
####################################################################################################################################
####################################################################################################################################


### --->>> Main function file that calculate the variables-of-interest for each country!! 
country.granger.test <- function(country.of.interest, input.data, num.lags = 1, num.COI = 10, num.impute = 5) {
  
  
  ### Main dataframe that will keep all the information
  final.df <- data.frame(chi.stat.AR = numeric(length = num.impute),
                         F.MEDIA =numeric(length = num.impute))
  
  ## Subset data for country-of-origin, and exclude total!
  test.data <- monthly.data[ monthly.data$GEO == country.of.interest & monthly.data$CITIZEN != "Total" , ]
  
  ## ----- DETERMINE WHICH COUNTRY-OF-ORIGINS ARE MOST SALIENT
  attach(test.data)
  salient.countries <- aggregate.data.frame(x = cbind(test.data$Total_positive, test.data$Rejected), 
                                            by = list(CITIZEN), FUN = sum, na.rm = TRUE)
  detach(test.data)
  colnames(salient.countries) <- c("CITIZEN", "Positive", "Rejected")
  salient.countries$Total <- salient.countries$Positive + salient.countries$Rejected
  salient.countries <- salient.countries[ order(salient.countries$Total, decreasing = TRUE) , ]
  
  ## Keep only those countries most salient
  test.data <- test.data[ test.data$CITIZEN %in% salient.countries$CITIZEN[1:num.COI] ,  ]
  test.data$CITIZEN <- factor(x = test.data$CITIZEN, levels = unique(test.data$CITIZEN))
  
  # Specify which columns are relevant
  keep.col <- c("TIME", "GEO", "CITIZEN", "Received", "Accept.Rate", "Tone.Rel", "NumArticles.Rel")
  test.data <- test.data[ , colnames(test.data) %in% keep.col ]
  
  ## Determine if we should proceed: is entire column full of NA values?
  if ( all(is.na(test.data$Accept.Rate)) ) {
    return(final.df)
  } else { #--- otherwise, proceed....
    
    # Generate multiple imputations of this dataset
    amelia.data <- amelia(test.data, # dataframe
                          m = num.impute, # number of imputations
                          ts = "TIME", # time-series element
                          cs = "CITIZEN", # cross-sectional element
                          idvars = c("GEO"), # identification variable 
                          lags = 1) # Number of lag-terms to include 
    
    ## Tone ~ Apps ???
    for (i in 1:num.impute) {
      ## Change time series unites
      use.data <- amelia.data$imputations[[i]]
      #time.vec <- seq(from = min(use.data$TIME),
      #                to = max(use.data$TIME),
      #                by = "month")
      time.vec <- sort(unique(use.data$TIME))
      time.cbind <- cbind(1:length(time.vec), time.vec)
      use.data$TIME = sapply(X = use.data$TIME, function(X) {time.cbind[ which(X == time.cbind[,2]) , 1]})
      
      ## Create panel data structure
      panel.data <- pdata.frame(x = use.data, 
                                index = c("CITIZEN", "TIME"), 
                                drop.index = TRUE, 
                                row.names = TRUE)
      ### Media --->>> Acceptance rate
      plm.model.1.AR <- plm(formula = Accept.Rate ~ -1 + plm::lag(x = Accept.Rate, k = 1:num.lags), 
                            data = panel.data, 
                            effect = "individual", model = "pooling")
      #summary(plm.model.1.AR)
      plm.model.2.AR <- plm(formula = Accept.Rate ~ -1 + plm::lag(x = Accept.Rate, k = 1:num.lags) + 
                              lag(x = Tone.Rel, k = 1:num.lags), 
                            data = panel.data,
                            effect = "individual", model = "pooling")
      #summary(plm.model.2.AR)
      wtest.AR <- waldtest(plm.model.1.AR, plm.model.2.AR)
      # Save the data
      final.df$chi.stat.AR[i] = min(wtest.AR$Chisq, na.rm = TRUE)
      
      ### Acceptance rate --->>> Media
      attach(amelia.data$imputations[[i]])
      agg.accept <- aggregate.data.frame(x = cbind(Accept.Rate, Tone.Rel),
                                         by = list(TIME), FUN = mean)
      detach(amelia.data$imputations[[i]])
      accept.ts <- ts(agg.accept$Accept.Rate[order(agg.accept$Group.1)], frequency = 12)
      media.ts <- ts(agg.accept$Tone.Rel[order(agg.accept$Group.1)], frequency = 12)
      plm.model.1.MEDIA <- dynlm(formula = media.ts ~ -1 + L(media.ts, 1:num.lags))
      #summary(plm.model.1.MEDIA)
      plm.model.2.MEDIA <- dynlm(formula = media.ts ~ -1 + L(media.ts, 1:num.lags) + L(accept.ts, 1:num.lags))
      #summary(plm.model.2.MEDIA)
      wtest.MEDIA <- waldtest(plm.model.1.MEDIA, plm.model.2.MEDIA)
      # Save the data
      final.df$F.MEDIA[i] = min(wtest.MEDIA$F, na.rm = TRUE)
      
    } ## End calculation of imputation
    
    
    return(final.df)
    
  } # End checking if column only has NA values
} ## End main-function calculating main variables for analysis

#---------------------------------------------------------------------------------------------

### --->>> Write a wrapper function that automatically computes results for all countries
wrapper.granger.result <- function(input.data, num.lags = 1, num.COI = 10, num.impute = 5) {
  
  # Define final results returned by function
  granger.results <- data.frame(country = unique(monthly.data$GEO),
                                AR.chi.NUM = numeric(length = length(unique(monthly.data$GEO))),
                                AR.chi.P = numeric(length = length(unique(monthly.data$GEO))),
                                MEDIA.F.NUM = numeric(length = length(unique(monthly.data$GEO))),
                                MEDIA.F.P = numeric(length = length(unique(monthly.data$GEO))))
  
  # Loop through every country
  for (i in 1:length(granger.results$country)) {
    # Calculate final results
    country.result <- country.granger.test(country.of.interest = granger.results$country[i], 
                                           input.data = monthly.data, num.lags = num.lags, 
                                           num.COI = num.COI, num.impute = num.impute)
    ## Check if country has valid acceptance rates!
    if ( sum(country.result$chi.stat.AR < 0.001) ) {
      granger.results$AR.chi.NUM[i] = NA
      granger.results$AR.chi.P[i] = NA
      granger.results$MEDIA.F.NUM[i] = NA
      granger.results$MEDIA.F.P[i] = NA
    } else {
      # Combine chi-square statistics
      val.AR <- micombine.chisquare(dk = country.result$chi.stat.AR, df = 1, display = FALSE)
      granger.results$AR.chi.NUM[i] = val.AR[1]
      granger.results$AR.chi.P[i] = val.AR[2]
      # Combine F statistics
      val.MEDIA <- micombine.F(Fvalues = country.result$F.MEDIA, df1 = 1, display = FALSE)
      granger.results$MEDIA.F.NUM[i] = val.MEDIA[1]
      granger.results$MEDIA.F.P[i] = val.MEDIA[2]
      
    } # End checking if country has values to impute!
  } # End looping through every country
  
  return(granger.results)
} ## End wrapper function for granger causality analysis

#---------------------------------------------------------------------------------------------

## "change.odds.ratio" Function
change.odds.ratio <- function(x,delta,b) {
  top <- (1+x+delta) * (1-x)
  bottom <- (1-x-delta) * (1+x)
  val <- (top / bottom)**b
  return(100*(val-1))
} # End "change.odds.ratio" function

#---------------------------------------------------------------------------------------------



####################################################################################################################################
####################################################################################################################################
#    PRODUCE TABLES FOR PAPER
####################################################################################################################################
####################################################################################################################################



## Small function calculating interesting model results
recover.R2 <- function(model.MEDIA, model.MEDIA.DIFF, model.RECEIVED, model.CONTROLS, model.FULL, model.FULL.DIFF) {
  
  return.df <- data.frame(Marginal.R2 = numeric(length = 6),
                          Conditional.R2 = numeric(length = 6))
  
  ## MEDIA
  r.squared.results <- as.data.frame(t(sapply(X = model.MEDIA, FUN =r.squaredGLMM)))
  return.df$Marginal.R2[1] = poolR2(input.R = r.squared.results[ , 1 ])
  return.df$Conditional.R2[1] = poolR2(input.R= r.squared.results[ , 2 ])
  
  ## MEDIA-DIFF
  r.squared.results <- as.data.frame(t(sapply(X = model.MEDIA.DIFF, FUN =r.squaredGLMM)))
  return.df$Marginal.R2[2] = poolR2(input.R = r.squared.results[ , 1 ])
  return.df$Conditional.R2[2] = poolR2(input.R= r.squared.results[ , 2 ])
  
  ## RECEIVED
  r.squared.results <- as.data.frame(t(sapply(X = model.RECEIVED, FUN =r.squaredGLMM)))
  return.df$Marginal.R2[3] = poolR2(input.R = r.squared.results[ , 1 ])
  return.df$Conditional.R2[3] = poolR2(input.R= r.squared.results[ , 2 ])
  
  ## CONTROLS
  r.squared.results <- as.data.frame(t(sapply(X = model.CONTROLS, FUN =r.squaredGLMM)))
  return.df$Marginal.R2[4] = poolR2(input.R = r.squared.results[ , 1 ])
  return.df$Conditional.R2[4] = poolR2(input.R= r.squared.results[ , 2 ])
  
  ## FULL
  r.squared.results <- as.data.frame(t(sapply(X = model.FULL, FUN =r.squaredGLMM)))
  return.df$Marginal.R2[5] = poolR2(input.R = r.squared.results[ , 1 ])
  return.df$Conditional.R2[5] = poolR2(input.R= r.squared.results[ , 2 ])
  
  ## MEDIA-DIFF
  r.squared.results <- as.data.frame(t(sapply(X = model.FULL.DIFF, FUN =r.squaredGLMM)))
  return.df$Marginal.R2[6] = poolR2(input.R = r.squared.results[ , 1 ])
  return.df$Conditional.R2[6] = poolR2(input.R= r.squared.results[ , 2 ])
  
  #
  return(return.df)
} # End printing results

#---------------------------------------------------------------------------------------------

### --- >>> Convert coefficient to latex format (& with p-value stars)
p.COEF <- function(num, p.val) {
  #... determine whether the number gets a star
  if (p.val < 0.001) {
    end.part <- "$^{***}$"
  } else if (p.val < 0.01) {
    end.part <- "$^{**}$"
  } else if (p.val < 0.05) {
    end.part <- "$^{*}$"
  } else {
    end.part <- ""
  } #... end determining whether number gets a star
  
  #... now print results
  new.num <- sub(pattern = "-", 
                 replacement = "$-$", 
                 x = as.character(round(x = num, digits = 3)))
  #... return results
  return(paste0(new.num, end.part))
} #... end function that converts number to latex

#---------------------------------------------------------------------------------------------



### --- >>> PRINT REGRESSION RESULTS FOR LATEX
print.Regressions <- function(with.intercept = FALSE,
                              regressions.MEDIA, regressions.MEDIA.DIFF, 
                              regressions.RECEIVED, regressions.CONTROLS, 
                              regressions.FULL, regressions.FULL.DIFF) {
  
  ## Shift all results by 1 if there is an intercept
  c <- ifelse(test = with.intercept, yes = 1, no = 0)
  
  
  # Print main summary values of coefficients!!!
  r.MEDIA <- testEstimates(model = regressions.MEDIA, var.comp = TRUE)$estimates
  r.MEDIA.DIFF <- testEstimates(model = regressions.MEDIA.DIFF, var.comp = TRUE)$estimates
  r.RECEIVED <- testEstimates(model = regressions.RECEIVED, var.comp = TRUE)$estimates
  r.CONTROLS <- testEstimates(model = regressions.CONTROLS, var.comp = TRUE)$estimates
  r.FULL <- testEstimates(model = regressions.FULL, var.comp = TRUE)$estimates
  r.FULL.DIFF <- testEstimates(model = regressions.FULL.DIFF, var.comp = TRUE)$estimates
  
  ## Build summaries
  s.MEDIA <- lapply(X = regressions.MEDIA, FUN = summary)
  s.MEDIA.DIFF <- lapply(X = regressions.MEDIA.DIFF, FUN = summary)
  s.RECEIVED <- lapply(X = regressions.RECEIVED, FUN = summary)
  s.CONTROLS <- lapply(X = regressions.CONTROLS, FUN = summary)
  s.FULL <- lapply(X = regressions.FULL, FUN = summary)
  s.FULL.DIFF <- lapply(X = regressions.FULL.DIFF, FUN = summary)
  
  x <- sapply(X = s.MEDIA, FUN = function(X) X$AICtab)
  
  #... now the printing procedure begins
  fileConn <- file("output.txt", "w")
  
  ##### BUILD STRING #####
  
  ##-------------------------------------------------------## ADD BREAKER
  write(x = " \\footnotesize \\textbf{Media measures: ($\\textbf{X}_i^t - \\overline{\\textbf{X}}_i^t$)} \\\\ ", file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## Sentiment, refugee $-$ baseline
  write(x = "\\hspace{5mm} \\footnotesize Sentiment, refugee $-$ baseline   &", file = fileConn, append = TRUE)
  write(x = paste0("  ", "&   ", #... MEDIA
                   "\\footnotesize ", p.COEF(num = r.MEDIA.DIFF[3+c,1], p.val = r.MEDIA.DIFF[3+c,5]), "   &   ",
                   "  ", " &   ", #... RECEIVED
                   "  ", " &   ", #... CONTROLS
                   "  ", " &   ", #... FULL
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[3+c,1], p.val = r.FULL.DIFF[3+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", #... <empty>
                   " ", "&   ", #... MEDIA
                   "\\tiny (", round(x = r.MEDIA.DIFF[3+c,2], digits = 3), ")   &   ",
                   "  ", "&   ", #... RECEIVED
                   "  ", "&   ", #... CONTROLS
                   "  ", "&   ", #... FULL
                   "\\tiny (", round(x = r.FULL.DIFF[3+c,2], digits = 3), ")  \\\\[0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## Sentiment, refugee
  write(x = "\\hspace{5mm} \\footnotesize Sentiment, refugee   &", file = fileConn, append = TRUE)
  write(x = paste0("\\footnotesize ", p.COEF(num = r.MEDIA[3+c,1], p.val = r.MEDIA[3+c,5]), " &   ",
                   "  ", "   &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[3+c,1], p.val = r.FULL[3+c,5]), " &   ",
                   "  ", "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "\\tiny (", round(x = r.MEDIA[3+c,2], digits = 3), ")   &   ",
                   " ", "&   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.FULL[3+c,2], digits = 3), ")   &   ",
                   "  \\\\[0.8ex] "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## Attention
  write(x = "\\hspace{5mm} \\footnotesize Attention   &", file = fileConn, append = TRUE)
  write(x = paste0("\\footnotesize ", p.COEF(num = r.MEDIA[1+c,1], p.val = r.MEDIA[1+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.MEDIA.DIFF[1+c,1], p.val = r.MEDIA.DIFF[1+c,5]), " &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[1+c,1], p.val = r.FULL[1+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[1+c,1], p.val = r.FULL.DIFF[1+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "\\tiny (", round(x = r.MEDIA[1+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.MEDIA.DIFF[1+c,2], digits = 3), ")   &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.FULL[1+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL.DIFF[1+c,2], digits = 3), ")   \\\\[0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## ADD BREAKER
  write(x = " \\footnotesize \\textbf{Media measures: $\\overline{\\textbf{X}}_i$} \\\\ ", file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## Sentiment, refugee  
  write(x = "\\hspace{5mm} \\footnotesize Sentiment, refugee   &", file = fileConn, append = TRUE)
  write(x = paste0("\\footnotesize ", p.COEF(num = r.MEDIA[4+c,1], p.val = r.MEDIA[4+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.MEDIA.DIFF[4+c,1], p.val = r.MEDIA.DIFF[4+c,5]), " &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[4+c,1], p.val = r.FULL[4+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[4+c,1], p.val = r.FULL.DIFF[4+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "\\tiny (", round(x = r.MEDIA[4+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.MEDIA.DIFF[4+c,2], digits = 3), ")   &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.FULL[4+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL.DIFF[4+c,2], digits = 3), ")   \\\\[0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## Attention  
  write(x = "\\hspace{5mm} \\footnotesize Attention   &", file = fileConn, append = TRUE)
  write(x = paste0("\\footnotesize ", p.COEF(num = r.MEDIA[2+c,1], p.val = r.MEDIA[2+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.MEDIA.DIFF[2+c,1], p.val = r.MEDIA.DIFF[2+c,5]), " &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[2+c,1], p.val = r.FULL[2+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[2+c,1], p.val = r.FULL.DIFF[2+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "\\tiny (", round(x = r.MEDIA[2+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.MEDIA.DIFF[2+c,2], digits = 3), ")   &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.FULL[2+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL.DIFF[2+c,2], digits = 3), ")   \\\\[0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## ADD BREAKER
  write(x = " \\footnotesize \\textbf{Time-dependent controls: ($\\textbf{X}_i^t - \\overline{\\textbf{X}}_i^t$)} \\\\ ", file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## Apps. received  
  write(x = "\\hspace{5mm} \\footnotesize Apps. received   &", file = fileConn, append = TRUE)
  write(x = paste0("  ", " &   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.RECEIVED[1+c,1], p.val = r.RECEIVED[1+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.CONTROLS[1+c,1], p.val = r.CONTROLS[1+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[5+c,1], p.val = r.FULL[5+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[5+c,1], p.val = r.FULL.DIFF[5+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.RECEIVED[1+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.CONTROLS[1+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL[5+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL.DIFF[5+c,2], digits = 3), ")   \\\\[0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## GDP  
  write(x = "\\hspace{5mm} \\footnotesize GDP   &", file = fileConn, append = TRUE)
  write(x = paste0("  ", " &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.CONTROLS[3+c,1], p.val = r.CONTROLS[3+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[7+c,1], p.val = r.FULL[7+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[7+c,1], p.val = r.FULL.DIFF[7+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "  ", "&   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.CONTROLS[3+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL[7+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL.DIFF[7+c,2], digits = 3), ")   \\\\[0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## Gov. debt  
  write(x = "\\hspace{5mm} \\footnotesize Gov. debt   &", file = fileConn, append = TRUE)
  write(x = paste0("  ", " &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.CONTROLS[5+c,1], p.val = r.CONTROLS[5+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[9+c,1], p.val = r.FULL[9+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[9+c,1], p.val = r.FULL.DIFF[9+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "  ", "&   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.CONTROLS[5+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL[9+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL.DIFF[9+c,2], digits = 3), ")   \\\\[0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## Unemployment rates
  write(x = "\\hspace{5mm} \\footnotesize Unemployment   &", file = fileConn, append = TRUE)
  write(x = paste0("  ", " &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.CONTROLS[7+c,1], p.val = r.CONTROLS[7+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[11+c,1], p.val = r.FULL[11+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[11+c,1], p.val = r.FULL.DIFF[11+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "  ", "&   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.CONTROLS[7+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL[11+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL.DIFF[11+c,2], digits = 3), ")   \\\\[0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## CPI
  write(x = "\\hspace{5mm} \\footnotesize CPI   &", file = fileConn, append = TRUE)
  write(x = paste0("  ", " &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.CONTROLS[9+c,1], p.val = r.CONTROLS[9+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[13+c,1], p.val = r.FULL[13+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[13+c,1], p.val = r.FULL.DIFF[13+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "  ", "&   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.CONTROLS[9+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL[13+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL.DIFF[13+c,2], digits = 3), ")   \\\\[0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## Gov. ideology (left < 0 < right) 
  write(x = "\\hspace{5mm} \\footnotesize Gov. ideology (left < 0 < right)   &", file = fileConn, append = TRUE)
  write(x = paste0("  ", " &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.CONTROLS[11+c,1], p.val = r.CONTROLS[11+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[15+c,1], p.val = r.FULL[15+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[15+c,1], p.val = r.FULL.DIFF[15+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "  ", "&   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.CONTROLS[11+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL[15+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL.DIFF[15+c,2], digits = 3), ")   \\\\[0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## Gov. ideology (left < 0 < right) 
  write(x = "\\hspace{5mm} \\footnotesize Press freedom index   &", file = fileConn, append = TRUE)
  write(x = paste0("  ", " &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.CONTROLS[13+c,1], p.val = r.CONTROLS[13+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[17+c,1], p.val = r.FULL[17+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[17+c,1], p.val = r.FULL.DIFF[17+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "  ", "&   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.CONTROLS[13+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL[17+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL.DIFF[17+c,2], digits = 3), ")   \\\\[0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## ADD BREAKER
  write(x = " \\footnotesize \\textbf{Time-invariant controls: $\\overline{\\textbf{X}}_i$} \\\\ ", file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## Apps. received  
  write(x = "\\hspace{5mm} \\footnotesize Apps. received   &", file = fileConn, append = TRUE)
  write(x = paste0("  ", " &   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.RECEIVED[2+c,1], p.val = r.RECEIVED[2+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.CONTROLS[2+c,1], p.val = r.CONTROLS[2+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[6+c,1], p.val = r.FULL[6+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[6+c,1], p.val = r.FULL.DIFF[6+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.RECEIVED[2+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.CONTROLS[2+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL[6+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL.DIFF[6+c,2], digits = 3), ")   \\\\[0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## GDP  
  write(x = "\\hspace{5mm} \\footnotesize GDP   &", file = fileConn, append = TRUE)
  write(x = paste0("  ", " &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.CONTROLS[4+c,1], p.val = r.CONTROLS[4+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[8+c,1], p.val = r.FULL[8+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[8+c,1], p.val = r.FULL.DIFF[8+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "  ", "&   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.CONTROLS[4+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL[8+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL.DIFF[8+c,2], digits = 3), ")   \\\\[0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## Gov. debt  
  write(x = "\\hspace{5mm} \\footnotesize Gov. debt   &", file = fileConn, append = TRUE)
  write(x = paste0("  ", " &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.CONTROLS[6+c,1], p.val = r.CONTROLS[6+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[10+c,1], p.val = r.FULL[10+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[10+c,1], p.val = r.FULL.DIFF[10+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "  ", "&   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.CONTROLS[6+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL[10+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL.DIFF[10+c,2], digits = 3), ")   \\\\[0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## Unemployment rates
  write(x = "\\hspace{5mm} \\footnotesize Unemployment   &", file = fileConn, append = TRUE)
  write(x = paste0("  ", " &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.CONTROLS[8+c,1], p.val = r.CONTROLS[8+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[12+c,1], p.val = r.FULL[12+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[12+c,1], p.val = r.FULL.DIFF[12+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "  ", "&   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.CONTROLS[8+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL[12+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL.DIFF[12+c,2], digits = 3), ")   \\\\[0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## CPI
  write(x = "\\hspace{5mm} \\footnotesize CPI   &", file = fileConn, append = TRUE)
  write(x = paste0("  ", " &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.CONTROLS[10+c,1], p.val = r.CONTROLS[10+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[14+c,1], p.val = r.FULL[14+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[14+c,1], p.val = r.FULL.DIFF[14+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "  ", "&   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.CONTROLS[10+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL[14+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL.DIFF[14+c,2], digits = 3), ")   \\\\   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## CPI
  write(x = "\\hspace{5mm} \\footnotesize Gov. ideology (left < 0 < right)   &", file = fileConn, append = TRUE)
  write(x = paste0("  ", " &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.CONTROLS[12+c,1], p.val = r.CONTROLS[12+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[16+c,1], p.val = r.FULL[16+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[16+c,1], p.val = r.FULL.DIFF[16+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "  ", "&   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.CONTROLS[12+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL[16+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL.DIFF[16+c,2], digits = 3), ")   \\\\   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## CPI
  write(x = "\\hspace{5mm} \\footnotesize Press freedom index   &", file = fileConn, append = TRUE)
  write(x = paste0("  ", " &   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\footnotesize ", p.COEF(num = r.CONTROLS[14+c,1], p.val = r.CONTROLS[14+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL[18+c,1], p.val = r.FULL[18+c,5]), " &   ",
                   "\\footnotesize ", p.COEF(num = r.FULL.DIFF[18+c,1], p.val = r.FULL.DIFF[18+c,5]), "  \\\\[-0.8ex]   "),
        file = fileConn, append = TRUE)
  write(x = paste0("  &  ", 
                   "  ", "&   ",
                   "  ", "&   ",
                   "  ", "&   ",
                   "\\tiny (", round(x = r.CONTROLS[14+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL[18+c,2], digits = 3), ")   &   ",
                   "\\tiny (", round(x = r.FULL.DIFF[18+c,2], digits = 3), ")   \\\\   "),
        file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  ##-------------------------------------------------------## ADD BREAKER
  write(x = "  " , file = fileConn, append = TRUE)
  write(x = " \\hline ", file = fileConn, append = TRUE)
  write(x = " \\hline \\\\[-1.8ex] ", file = fileConn, append = TRUE)
  write(x = "  " , file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  n <- regressions.MEDIA$imp1@devcomp$dims[[1]]
  
  ##-------------------------------------------------------## Num. observations
  write(x = "\\hspace{5mm} \\footnotesize Num. observations   &  ", file = fileConn, append = TRUE)
  write(x = paste0("\\footnotesize ", n, " &   ",
                   "\\footnotesize ", n, " &   ",
                   "\\footnotesize ", n, " &   ",
                   "\\footnotesize ", n, " &   ",
                   "\\footnotesize ", n, " &   ",
                   "\\footnotesize ", n, "   \\\\  "),
        file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  #..
  #..
  #..
  
  r <- recover.R2(model.MEDIA = regressions.MEDIA, model.MEDIA.DIFF = regressions.MEDIA.DIFF, 
                  model.RECEIVED = regressions.RECEIVED, model.CONTROLS = regressions.CONTROLS, 
                  model.FULL = regressions.FULL, model.FULL.DIFF = regressions.FULL.DIFF)
  
  ##-------------------------------------------------------## Marginal $R^2$
  write(x = "\\hspace{5mm} \\footnotesize Marginal $R^2$   &  ", file = fileConn, append = TRUE)
  write(x = paste0("\\footnotesize ", round(x = r$Marginal.R2[1], digits = 2), " &   ",
                   "\\footnotesize ", round(x = r$Marginal.R2[2], digits = 2), " &   ",
                   "\\footnotesize ", round(x = r$Marginal.R2[3], digits = 2), " &   ",
                   "\\footnotesize ", round(x = r$Marginal.R2[4], digits = 2), " &   ",
                   "\\footnotesize ", round(x = r$Marginal.R2[5], digits = 2), " &   ",
                   "\\footnotesize ", round(x = r$Marginal.R2[6], digits = 2),"   \\\\  "),
        file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  ##-------------------------------------------------------## Conditional $R^2$
  write(x = "\\hspace{5mm} \\footnotesize Conditional $R^2$   &  ", file = fileConn, append = TRUE)
  write(x = paste0("\\footnotesize ", round(x = r$Conditional.R2[1], digits = 2), " &   ",
                   "\\footnotesize ", round(x = r$Conditional.R2[2], digits = 2), " &   ",
                   "\\footnotesize ", round(x = r$Conditional.R2[3], digits = 2), " &   ",
                   "\\footnotesize ", round(x = r$Conditional.R2[4], digits = 2), " &   ",
                   "\\footnotesize ", round(x = r$Conditional.R2[5], digits = 2), " &   ",
                   "\\footnotesize ", round(x = r$Conditional.R2[6], digits = 2),"   \\\\  "),
        file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  ##-------------------------------------------------------## AICc
  write(x = "\\hspace{5mm} \\footnotesize AICc   &  ", file = fileConn, append = TRUE)
  write(x = paste0("\\footnotesize ", round(x = mean(sapply(X = s.MEDIA, FUN = function(X) X$AICtab)[1,]), digits = 0), " &   ",
                   "\\footnotesize ", round(x = mean(sapply(X = s.MEDIA.DIFF, FUN = function(X) X$AICtab)[1,]), digits = 0), " &   ",
                   "\\footnotesize ", round(x = mean(sapply(X = s.RECEIVED, FUN = function(X) X$AICtab)[1,]), digits = 0), " &   ",
                   "\\footnotesize ", round(x = mean(sapply(X = s.CONTROLS, FUN = function(X) X$AICtab)[1,]), digits = 0), " &   ",
                   "\\footnotesize ", round(x = mean(sapply(X = s.FULL, FUN = function(X) X$AICtab)[1,]), digits = 0), " &   ",
                   "\\footnotesize ", round(x = mean(sapply(X = s.FULL.DIFF, FUN = function(X) X$AICtab)[1,]), digits = 0),"   \\\\  "),
        file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  ##-------------------------------------------------------## BIC
  write(x = "\\hspace{5mm} \\footnotesize BIC   &", file = fileConn, append = TRUE)
  write(x = paste0("\\footnotesize ", round(x = mean(sapply(X = s.MEDIA, FUN = function(X) X$AICtab)[2,]), digits = 0), " &   ",
                   "\\footnotesize ", round(x = mean(sapply(X = s.MEDIA.DIFF, FUN = function(X) X$AICtab)[2,]), digits = 0), " &   ",
                   "\\footnotesize ", round(x = mean(sapply(X = s.RECEIVED, FUN = function(X) X$AICtab)[2,]), digits = 0), " &   ",
                   "\\footnotesize ", round(x = mean(sapply(X = s.CONTROLS, FUN = function(X) X$AICtab)[2,]), digits = 0), " &   ",
                   "\\footnotesize ", round(x = mean(sapply(X = s.FULL, FUN = function(X) X$AICtab)[2,]), digits = 0), " &   ",
                   "\\footnotesize ", round(x = mean(sapply(X = s.FULL.DIFF, FUN = function(X) X$AICtab)[2,]), digits = 0),"   \\\\  "),
        file = fileConn, append = TRUE)
  ##-------------------------------------------------------##
  
  
  
  ### End by closing the text file
  close(fileConn)
  
}#... end printing latex results for latex



#---------------------------------------------------------------------------------------------


### Function for producing data about press freedom
# Data comes from here: https://rsf.org/en/ranking
produce.press.data <- function() {

  Austria.press <- rev(c(13.47,#	(011)
                         13.18,#	(007)
                         10.85,#	(012)
                         10.01,#	(012)
                         9.40	,#(005)
                         -8.00,#	(007)
                         0.50,#	(013)
                         3.00,#	(014)
                         3.50	,#(016)
                         4.25,#	(016)
                         4.50,#	(016)
                         2.50,#	(017)
                         3.25,#	(016)
                         2.75,#	(026)
                         7.50))
  
  Belgium.press <- rev(c(12.75,#	(013)
                         14.18,#	(015)
                         11.98,#	(023)
                         12.80,#	(021)
                         12.94,#	(020)
                         -2.00,#	(014)
                         4.00,#	(011)
                         2.50,#	(007)
                         3.00,#	(005)
                         1.50,#	(014)
                         4.00,#	(018)
                         4.00,#	(022)
                         4.00,#	(007)
                         1.17,#	(012)
                         3.50))
  
  Czech.press <- rev(c(16.91,#	(021)
                       16.66,#	(013)
                       11.62,#	(013)
                       10.07,#	(016)
                       10.17,#	(014)
                       -5.00,#	(023)
                       7.50,#	(024)
                       5.00,#	(016)
                       4.00,#	(014)
                       4.00,#	(005)
                       0.75,#	(009)
                       1.00,#	(019)
                       3.50,#	(012)
                       2.50,#	(041)
                       11.25))
  
  Denmark.press <- rev(c(10.36,#	(004)
                         8.89,#	(003)
                         8.24,#	(007)
                         7.43,#	(006)
                         7.08,#	(010)
                         -5.67,#	(011)
                         2.50,#	(001)
                         0.00,#	(014)
                         3.50,#	(008)
                         2.00,#	(019)
                         5.00,#	(001)
                         0.50,#	(001)
                         0.50,#	(005)
                         1.00,#	(010)
                         3.00))
  
  Estonia.press <- rev(c(13.55,#	(014)
                         14.31,#	(010)
                         11.19,#	(011)
                         9.63,#	(011)
                         9.26,#	(003)
                         -9.00,#	(009)
                         2.00,#	(006)
                         0.50,#	(004)
                         2.00,#	(003)
                         1.00,#	(006)
                         2.00,#	(011)
                         1.50,#	(011)
                         2.00,#	(012)
                         2.50,
                         2.5))
  
  Finland.press <- rev(c(8.92,#	(001)
                         8.59,#	(001)
                         7.52,#	(001)
                         6.40,#	(001)
                         6.38,#	(001)
                         -10.00,#	(001)
                         0.00,#	(001)
                         0.00,#	(004)
                         2.00,#	(005)
                         1.50,#	(001)
                         0.50,#	(001)
                         0.50,#	(001)
                         0.50,#	(001)
                         0.50,#	(001)
                         0.50))
  
  France.press <- rev(c(22.24,#	(045)
                        23.83,#	(038)
                        21.15,#	(039)
                        21.89,#	(037)
                        21.60,#	(038)
                        9.50,#	(044)
                        13.38,#	(043)
                        10.67,#	(035)
                        7.67,#	(031)
                        9.75,#	(035)
                        9.00,#	(030)
                        6.25,#	(019)
                        3.50,#	(026)
                        4.17,#	(011)
                        3.25))
  
  Germany.press <- rev(c(14.97,#	(016)
                         14.80,#	(012)
                         11.47,#	(014)
                         10.23,#	(017)
                         10.24,#	(016)
                         -3.00,#	(017)
                         4.25,#	(018)
                         3.50,#	(020)
                         4.50,#	(020)
                         5.75,#	(023)
                         5.50,#	(018)
                         4.00,#	(011)
                         2.00,#	(008)
                         1.33,#	(007)
                         1.50))
  
  
  Greece.press <- rev(c(30.89,#	(089)
                        30.35,#	(091)
                        31.01,#	(099)
                        31.33,#	(084)
                        28.46,#	(070)
                        24.00,#	(070)
                        19.00,#	(035)
                        9.00,#	(031)
                        7.50,#	(030)
                        9.25,#	(032)
                        8.00,#	(018)
                        4.00,#	(033)
                        7.00,#	(031)
                        6.00,#	(019)
                        5.00))
  
  Hungary.press <- rev(c(29.01,#	(067)
                         28.17,#	(065)
                         27.44,#	(064)
                         26.73,#	(056)
                         26.09,#	(040)
                         10.00,#	(023)
                         7.50,#	(025)
                         5.50,#	(023)
                         5.50,#	(017)
                         4.50,#	(010)
                         3.00,#	(012)
                         2.00,#	(028)
                         6.00,#	(021)
                         3.33,#	(024)
                         6.50))
  Italy.press <- rev(c(26.26,#	(077)
                       28.93,#	(073)
                       27.94,#	(049)
                       23.75,#	(057)
                       26.11,#	(061)
                       19.67,#	(049)
                       15.00,#	(049)
                       12.14,#	(044)
                       8.42	,#(035)
                       11.25,#	(040)
                       9.90	,#(042)
                       8.67,#	(039)
                       9.00,#	(053)
                       9.75,#	(040)
                       11.00))
  Latvia.press <- rev(c(18.62,#	(024)
                        17.38,#	(028)
                        18.12,#	(037)
                        21.10,#	(039)
                        22.89,#	(050)
                        15.00,#	(030)
                        8.50,#	(013)
                        3.00,#	(007)
                        3.00,#	(012)
                        3.50,#	(010)
                        3.00,#	(016)
                        2.50,#	(010)
                        1.00,#	(011)
                        2.25,
                        NA))
  
  
  Lithuania.press <- rev(c(21.37,#	(035)
                           19.95,#	(031)
                           18.80,#	(032)
                           19.20,#
                           18.24,#
                           4.00,#	(011)
                           2.50,#	(010)
                           2.25,#	(016)
                           4.00,#	(023)
                           7.00,#	(027)
                           6.50,#	(021)
                           4.50,#	(016)
                           3.00,#	(017)
                           2.83,
                           NA))
  
  Netherlands.press <- rev(c(11.28,
                             8.76,#	(004)
                             9.22,#	(002)
                             6.46,#	(002)
                             6.48,#	(003)
                             -9.00,#	(001)
                             0.00,#	(007)
                             1.00,#	(016)
                             4.00,#	(012)
                             3.50,#	(001)
                             0.50,#	(001)
                             0.50,#	(001)
                             0.50,#	(001)
                             0.50,#	(001)
                             0.50))
  
  Norway.press <- rev(c(7.6,#	(003)
                        8.79,#	(002)
                        7.75,#	(003)
                        6.52,#	(003)
                        6.52,#	(001)
                        -10.00,#	(001)
                        0.00,#	(001)
                        0.00,#	(001)
                        1.50,#	(001)
                        0.75,#	(006)
                        2.00,#	(001)
                        0.50,#	(001)
                        0.50,#	(001)
                        0.50,#	(001)
                        0.50))
  Poland.press <- rev(c(26.47,#	(047)
                        23.89,#	(018)
                        12.71,#	(019)
                        11.03,#	(022)
                        13.11,#	(024)
                        -0.67,#	(032)
                        8.88,#	(037)
                        9.50,#	(047)
                        9.00,#	(056)
                        18.50,#	(058)
                        14.00,#	(053)
                        12.50,#	(032)
                        6.83,#	(033)
                        6.17,#	(029)
                        7.75))
  Portugal.press <- rev(c(15.77,#	(023)
                          17.27,#	(026)
                          17.11,#	(030)
                          17.73,#	(028)
                          16.75,#	(033)
                          5.33,#	(040)
                          12.36,#	(030)
                          8.00,#	(016)
                          4.00,#	(008)
                          2.00,#	(010)
                          3.00,#	(023)
                          4.83,#	(025)
                          4.50,#	(028)
                          5.17,#	(007)
                          1.50))
  Spain.press <- rev(c(18.69,#	(034)
                       19.92,#	(033)
                       19.95,#	(035)
                       20.63,#	(036)
                       20.50,#	(039)
                       9.75	,#(039)
                       12.25,#	(044)
                       11.00,#	(036)
                       8.00	,#(033)
                       10.25,#	(041)
                       10.00,#	(040)
                       8.33	,#(039)
                       9.00	,#(042)
                       7.67	,#(029)
                       7.75))
  Sweden.press <- rev(c(8.27,#	(008)
                        12.33,#	(005)
                        9.47,#	(010)
                        8.98,#	(010)
                        9.23,#	(012)
                        -5.50,#	(001)
                        0.00,#	(001)
                        0.00,#	(007)
                        3.00,#	(005)
                        1.50,#	(014)
                        4.00,#	(012)
                        2.00,#	(011)
                        2.00,#	(009)
                        1.50,#	(007)
                        1.50))
  UK.press <- rev(c(22.26,#	(038)
                    21.70,#	(034)
                    20.00,#	(033)
                    19.93,#	(029)
                    16.89,#	(028)
                    2.00,#	(019)
                    6.00,#	(020)
                    4.00,#	(023)
                    5.50,#	(024)
                    8.25,#	(027)
                    6.50,#	(024)
                    5.17,#	(028)
                    6.00,#	(027)
                    4.25,#	(021)
                    6.00))
  
  ### European countries
  European.countries <- c("Austria", "Belgium","Czech Republic","Denmark","Estonia","Finland",
                          "France","Germany","Greece","Hungary","Italy","Latvia","Lithuania",
                          "Netherlands","Norway","Poland","Portugal","Spain","Sweden","United Kingdom")
  
  ### Year-quarter
  time.vec <- as.yearqtr(c("2002 Q1","2002 Q2","2002 Q3", "2002 Q4", "2003 Q1", "2003 Q2", "2003 Q3", "2003 Q4", "2004 Q1", "2004 Q2", "2004 Q3", "2004 Q4", "2005 Q1", "2005 Q2", "2005 Q3",
                "2005 Q4", "2006 Q1", "2006 Q2", "2006 Q3", "2006 Q4", "2007 Q1", "2007 Q2", "2007 Q3", "2007 Q4", "2008 Q1", "2008 Q2", "2008 Q3", "2008 Q4", "2009 Q1", "2009 Q2",
                "2009 Q3", "2009 Q4", "2010 Q1", "2010 Q2", "2010 Q3", "2010 Q4", "2011 Q1", "2011 Q2", "2011 Q3", "2011 Q4", "2012 Q1", "2012 Q2", "2012 Q3", "2012 Q4", "2013 Q1",
                "2013 Q2", "2013 Q3", "2013 Q4", "2014 Q1", "2014 Q2", "2014 Q3", "2014 Q4", "2015 Q1", "2015 Q2", "2015 Q3", "2015 Q4", "2016 Q1", "2016 Q2", "2016 Q3", "2016 Q4"))
  
  ## Combine Qarters and years
  return.data <- data.frame(GEO = rep(x = European.countries, each = length(time.vec)),
                           TIME = rep(x = time.vec, times = length(European.countries)))
  ## Add empty vector... for adding press values later
  return.data$Press.Freedom <- numeric(length = length(return.data$GEO))
  
  ## Loop through every European country 
  return.data[ return.data$GEO == European.countries[1]  , "Press.Freedom"] = rep(x = Austria.press, each = 4) #... one for each quarter
  return.data[ return.data$GEO == European.countries[2]  , "Press.Freedom"] = rep(x = Belgium.press, each = 4) 
  return.data[ return.data$GEO == European.countries[3]  , "Press.Freedom"] = rep(x = Czech.press, each = 4) 
  return.data[ return.data$GEO == European.countries[4]  , "Press.Freedom"] = rep(x = Denmark.press, each = 4) 
  return.data[ return.data$GEO == European.countries[5]  , "Press.Freedom"] = rep(x = Estonia.press, each = 4) 
  return.data[ return.data$GEO == European.countries[6]  , "Press.Freedom"] = rep(x = Finland.press, each = 4) 
  return.data[ return.data$GEO == European.countries[7]  , "Press.Freedom"] = rep(x = France.press, each = 4) 
  return.data[ return.data$GEO == European.countries[8]  , "Press.Freedom"] = rep(x = Germany.press, each = 4) 
  return.data[ return.data$GEO == European.countries[9]  , "Press.Freedom"] = rep(x = Greece.press, each = 4) 
  return.data[ return.data$GEO == European.countries[10]  , "Press.Freedom"] = rep(x = Hungary.press, each = 4) 
  return.data[ return.data$GEO == European.countries[11]  , "Press.Freedom"] = rep(x = Italy.press, each = 4) 
  return.data[ return.data$GEO == European.countries[12]  , "Press.Freedom"] = rep(x = Latvia.press, each = 4) 
  return.data[ return.data$GEO == European.countries[13]  , "Press.Freedom"] = rep(x = Lithuania.press, each = 4) 
  return.data[ return.data$GEO == European.countries[14]  , "Press.Freedom"] = rep(x = Netherlands.press, each = 4) 
  return.data[ return.data$GEO == European.countries[15]  , "Press.Freedom"] = rep(x = Norway.press, each = 4) 
  return.data[ return.data$GEO == European.countries[16]  , "Press.Freedom"] = rep(x = Poland.press, each = 4) 
  return.data[ return.data$GEO == European.countries[17]  , "Press.Freedom"] = rep(x = Portugal.press, each = 4) 
  return.data[ return.data$GEO == European.countries[18]  , "Press.Freedom"] = rep(x = Spain.press, each = 4) 
  return.data[ return.data$GEO == European.countries[19]  , "Press.Freedom"] = rep(x = Sweden.press, each = 4) 
  return.data[ return.data$GEO == European.countries[20]  , "Press.Freedom"] = rep(x = UK.press, each = 4) 
  
  return(return.data)
  
} #... end producing data about the freedom of press



#---------------------------------------------------------------------------------------------


## Add ideology scores to the dataset
# Data comes from here: https://manifesto-project.wzb.eu/
produce.Ideology.data <- function(Ideology) {
  
  ## Load the data
  Ideology$rel.Ideology <- (Ideology$absseat / Ideology$totseats) * Ideology$rile
  attach(Ideology)
  agg.Ideology <- aggregate.data.frame(x = rel.Ideology, by = list(countryname, date), FUN = sum, na.rm = TRUE)
  detach(Ideology)
  agg.Ideology$Year <- sapply(X = agg.Ideology$Group.2, FUN = function(X) substr(x = X, start = 1, stop = 4))
  agg.Ideology$Month <- sapply(X = agg.Ideology$Group.2, FUN = function(X) substr(x = X, start = 5, stop = 6))
  agg.Ideology$Quarter <- ceiling(as.numeric(agg.Ideology$Month) / 3)
  agg.Ideology$TIME <- apply(X = cbind(agg.Ideology$Year, agg.Ideology$Quarter), MARGIN = 1, FUN = function(X) as.yearqtr(x = paste0(X[1], " Q", X[2]), format = "%Y Q%q") )
  
  ## Clean the dataset
  agg.Ideology <- agg.Ideology[ agg.Ideology$TIME >= as.yearqtr("1995 Q1") , ]
  agg.Ideology <- agg.Ideology[ , c("Group.1", "x", "TIME")]
  colnames(agg.Ideology) = c("GEO", "Ideology", "TIME")
  agg.Ideology$TIME <- as.yearqtr(agg.Ideology$TIME)
  
  ## Now fill in the details
  
  ### European countries
  European.countries <- c("Austria", "Belgium","Czech Republic","Denmark","Estonia","Finland",
                          "France","Germany","Greece","Hungary","Italy","Latvia","Lithuania",
                          "Netherlands","Norway","Poland","Portugal","Spain","Sweden","United Kingdom")
  
  ### Year-quarter
  time.vec <- as.yearqtr(c("2000 Q1","2000 Q2","2000 Q3", "2000 Q4", 
                           "2001 Q1","2001 Q2","2001 Q3", "2001 Q4", 
                           "2002 Q1","2002 Q2","2002 Q3", "2002 Q4", 
                           "2003 Q1", "2003 Q2", "2003 Q3", "2003 Q4", 
                           "2004 Q1", "2004 Q2", "2004 Q3", "2004 Q4", 
                           "2005 Q1", "2005 Q2", "2005 Q3", "2005 Q4", 
                           "2006 Q1", "2006 Q2", "2006 Q3", "2006 Q4", 
                           "2007 Q1", "2007 Q2", "2007 Q3", "2007 Q4", 
                           "2008 Q1", "2008 Q2", "2008 Q3", "2008 Q4", 
                           "2009 Q1", "2009 Q2", "2009 Q3", "2009 Q4", 
                           "2010 Q1", "2010 Q2", "2010 Q3", "2010 Q4", 
                           "2011 Q1", "2011 Q2", "2011 Q3", "2011 Q4", 
                           "2012 Q1", "2012 Q2", "2012 Q3", "2012 Q4", 
                           "2013 Q1", "2013 Q2", "2013 Q3", "2013 Q4", 
                           "2014 Q1", "2014 Q2", "2014 Q3", "2014 Q4", 
                           "2015 Q1", "2015 Q2", "2015 Q3", "2015 Q4", 
                           "2016 Q1", "2016 Q2", "2016 Q3", "2016 Q4"))
  
  ## Combine Qarters and years
  Ideology.final.data <- data.frame(GEO = rep(x = European.countries, each = length(time.vec)),
                                    TIME = rep(x = time.vec, times = length(European.countries)))
  ## Add empty vector... for adding press values later
  Ideology.final.data$Ideology <- numeric(length = length(Ideology.final.data$GEO))
  Ideology.final.data$TIME <- as.yearqtr(Ideology.final.data$TIME)
  
  
  ## Try to build this into results
  for (i in 1:length(Ideology.final.data$GEO)) {
    # Extract the year.qrt
    yQ <- Ideology.final.data$TIME[i]
    # Extract country
    country <- as.character(Ideology.final.data$GEO[i])
    # Extract the unique set of election quarters for this country
    election.quarters <- unique(agg.Ideology[ agg.Ideology$GEO == country , "TIME"])
    # If the present time is in agg.Ideology, then add it!
    if (yQ %in% election.quarters) {
      Ideology.final.data$Ideology[i] = mean(agg.Ideology[ agg.Ideology$GEO == country & agg.Ideology$TIME == yQ , "Ideology" ])
    } else { #... otherwise, look for past elections to match the result
      #.. find date that is closest to it
      x <- as.numeric(election.quarters) - as.numeric(yQ)
      index.date <- which.max(x[x < 0])
      Ideology.final.data$Ideology[i] = mean(agg.Ideology[ agg.Ideology$GEO == country & agg.Ideology$TIME == election.quarters[index.date] , "Ideology" ])
    } #... end checking if this is or is not an election year
  } #... end looping through the entire data set to add Ideology scores
  
  ## Return the final dataset
  return(Ideology.final.data)
} ##... end adding the ideology scores to each dataset


#---------------------------------------------------------------------------------------------

## Stack values in order to produce the time-series plot for Figure 1(b)
generate.stack.values <- function(countries.keep, month.val, country.data) {
  return.vec <- c()
  for (c in countries.keep) {
    return.vec <- c(return.vec, country.data[country.data$GEO == c &
                                              country.data$TIME == month.val , ]$Total_positive)
  }
  return.vec <- c( return.vec , sum( na.omit(country.data[country.data$TIME == month.val , ]$Total_positive)) - 
                     sum(na.omit(return.vec)) )
  ## Now return the vector
  return.vec[is.na(return.vec)] = 0
  return(return.vec)
} ## End function generating data



#---------------------------------------------------------------------------------------------

## Function for producing asylum plots in the main text
produce.asylum.plot <- function(a, country, r.data) {
  ## Extract asylum data
  Asylum.data <- r.data[ r.data$CITIZEN == a , ]
  
 # pdf(file = paste0("Fig_2_", country, ".pdf"), width = 7, height = 4)

  ## Germany
  x <- Asylum.data[Asylum.data$GEO == country[1] , "TIME"]
  x <- x[!is.na(x)]
  y <- Asylum.data[Asylum.data$GEO == country[1], "Acceptance.Rate"]
  y <- y[!is.na(y)]
  if (length(x) > length(y)) {
    x <- x[-1]
  } else if (length(x) < length(y)) {
    y <- y[-1]
  }
  plot(x,y, ylim = c(0,1), axes = FALSE, ylab = "Asylum acceptance rate (per quarter)", xlab = "Time", main = "")
  lines(x,y,lwd=2)
  
  ## UK
  x <- Asylum.data[Asylum.data$GEO == country[2] , "TIME"]
  x <- x[!is.na(x)]
  y <- Asylum.data[Asylum.data$GEO == country[2], "Acceptance.Rate"]
  y <- y[!is.na(y)]
  if (length(x) > length(y)) {
    x <- x[-1]
  } else if (length(x) < length(y)) {
    y <- y[-1]
  }
  points(x,y,col = 'black', pch = 15)
  lines(x,y,col = 'black', lty = 3, lwd=2)
  
  ## France
  x <- Asylum.data[Asylum.data$GEO == country[3] , "TIME"]
  x <- x[!is.na(x)]
  y <- Asylum.data[Asylum.data$GEO == country[3], "Acceptance.Rate"]
  y <- y[!is.na(y)]
  if (length(x) > length(y)) {
    x <- x[-1]
  } else if (length(x) < length(y)) {
    y <- y[-1]
  }
  points(x,y, col = 'black', pch = 4)
  lines(x,y,col = 'black', lty = 4, lwd=2)
  
  
  ## Add everything ex post
  time <- unique(Asylum.data$TIME)
  time <- time[grep(pattern = "Q1", x = time)]
  axis(side = 2, at = c(0.0, 0.2, 0.4, 0.6, 0.8, 1.0), labels = c(0.0, 0.2, 0.4, 0.6, 0.8, 1.0))
  axis(side = 1, at = time, labels = 2002:2016)
  title(main = a)
  legend("bottomright", legend = country, col = "black", lty = c(1,3,4), pch = c(1,15,4), bty = "n")
  
  # dev.off()
  
} ## End "produce.asylum.plot" function

